 
#include "glodef.h"
#include "errlog.h" 
#include "grid3d.h" 

#include <stdlib.h>
#include <stdio.h>
#include <string.h>


/****************************************/
/* print_nodes */
/****************************************/

void print_nodes(grid3d_t *grid, const node_t *nds, char *fname)
{
    FILE *F;
    int i;
    const node_t *n;
    double x, y, z;

    F=fopen(fname, "wt");
    if(F==NULL)
    {
        perror("print_nodes");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    if((grid->NN==0)||(grid->nodes==NULL))
        goto close;

    if(nds==NULL)
        nds=grid->nodes;

    fprintf(F, "%6s %15s %15s %15s %6s \n", "No", "x", "y", "z", "bcond");
    for(i=0, n=nds; i<grid->NN; i++, n++)
    {
        x=n->coords[0];
        y=n->coords[1];
        z=n->coords[2];
        fprintf(F, "%6d %15.12g %15.12g %15.12g %6d\n", i, x, y, z, n->bcond);
    }

close:
    fclose(F);
    return;
}

/****************************************/
/* print_edges */

/****************************************/

void print_edges(grid3d_t *grid, char *fname)
{
    FILE *F;
    int i, j, max_enf;
    edge_t *e;
    double t[3];
    double l, T[3];

    F=fopen(fname, "wt");
    if(F==NULL)
    {
        perror("print_edges");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    if((grid->NE==0)||(grid->edges==NULL))
        goto close;

    fprintf(F, "%6s %6s %6s %6s %10s %10s %10s %10s %6s %6s\n",
        "No", "n0", "n1", "bcond", "length",
        "Tx", "Ty", "Tz", "nf", "faces");
    for(max_enf=i=0, e=grid->edges; i<grid->NE; i++, e++)
        max_enf=max(max_enf, e->nf);
    for(i=0, e=grid->edges; i<grid->NE; i++, e++)
    {
        l=edge_length(grid, e, t);
        for(j=0; j<3; j++)
            T[j]=t[j];
        fprintf(F, "%6d %6d %6d %6d %10.5g %10.5g %10.5g %10.5g %6d",
            i, e->nodes[0], e->nodes[1], e->bcond, l,
            T[0], T[1], T[2], e->nf);
        for(j=0; j<e->nf; j++)
            fprintf(F, " %6d", e->faces[j]);
        for(; j<max_enf; j++)
            fprintf(F, " %6d", -1);
        fputc('\n', F);
    }

close:
    fclose(F);
    return;
}

/****************************************/
/* print_faces */

/****************************************/

void print_faces(grid3d_t *grid, char *fname)
{
    FILE *F;
    int i;
    face_t *f;
    double n[3];
    double A, N[3];

    F=fopen(fname, "wt");
    if(F==NULL)
    {
        perror("print_faces");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    if((grid->NF==0)||(grid->faces==NULL))
        goto close;

    fprintf(F, "%6s %6s %6s %6s %6s %6s %6s %6s %6s %6s %10s %10s %10s %10s\n",
        "No", "n0", "n1", "n2", "e0", "e1", "e2", "T0", "T1", "bcond",
        "Area", "Nx", "Ny", "Nz");
    for(i=0, f=grid->faces; i<grid->NF; i++, f++)
    {
        A=face_area(grid, f, n);
        N[0]=n[0];
        N[1]=n[1];
        N[2]=n[2];
        fprintf(F, "%6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %10.5g %10.5g"
            "%10.5g %10.5g\n", i, f->nodes[0], f->nodes[1], f->nodes[2],
            f->edges[0], f->edges[1], f->edges[2], f->elems[0], f->elems[1],
            f->bcond, A, N[0], N[1], N[2]);
    }

close:
    fclose(F);
    return;
}

/****************************************/
/* print_elems */

/****************************************/

void print_elems(grid3d_t *grid, char *fname)
{
    FILE *F;
    int i;
    elem_t *T;
    double V;

    F=fopen(fname, "wt");
    if(F==NULL)
    {
        perror("print_elems");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    if((grid->NEL==0)||(grid->elems==NULL))
        goto close;

    fprintf(F, "%6s %6s %6s %6s %6s %6s %6s %6s %6s %10s\n",          \
            "No", "nodes", "", "", "", "faces", "", "", "", "Volume");
    for(i=0, T=grid->elems; i<grid->NEL; i++, T++)
    {
        V=elem_volume(grid, T);
        fprintf(F, "%6d %6d %6d %6d %6d %6d %6d %6d %6d %10.5g\n",          \
                i, T->nodes[0], T->nodes[1], T->nodes[2], T->nodes[3],          \
                T->faces[0], T->faces[1], T->faces[2], T->faces[3], V);
    }

close:
    fclose(F);
    return;
}

/****************************************************************************/
/*            SaveGrid and SaveGridBin */
/****************************************************************************/
/* NOTE: SaveGrid overwrites fname if it exists, while SaveGridBin appends */

/*  binary data at the end of fname or creates new file */

void SaveGrid(const grid3d_t *grid, const char *fname)
{
    FILE *F;
    int i;

    F=fopen(fname, "wt");
    if(F==NULL)
    {
        perror("SaveGrid");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    fprintf(F, "%d # number of nodes\n", grid->NN);
    for(i=0; i<grid->NN; i++)
    {
        double *p=grid->nodes[i].coords;
        fprintf(F, "%.19g %.19g %.19g\n",
            (double)p[0], (double)p[1], (double)p[2]);
    }
    fprintf(F, "\n");
    for(i=0; i<grid->NN; i++)
    {
        fprintf(F, "%u\n", grid->nodes[i].on_bdry);
    }
    fprintf(F, "\n");
    fprintf(F, "%d # number of elements\n", grid->NEL);
    for(i=0; i<grid->NEL; i++)
    {
        int *en=grid->elems[i].nodes;
        fprintf(F, "%d %d %d %d\n",
            en[0], en[1], en[2], en[3]);
    }
    fclose(F);
}

void SaveGridBin(const grid3d_t *grid, const char *fname)
{
    size_t ignore;
    FILE *F;

    if(grid==NULL)
        return;
    F=fopen(fname, "ab");
    if(F==NULL)
    {
        perror("SaveGridBin");
        error_msg("%s - cannot open file for writing.\n", fname);
    }
    fseek(F, 0, SEEK_END);
    if(ferror(F))
    {
        perror("SaveGridBin");
        error_msg("SaveGridBin: fseek failed");
    }
    ignore=fwrite(&grid->NN, sizeof(int), 1, F);
    {
        int i;
        double p[3];
        for(i=0; i<grid->NN; i++)
        {
            p[0]=grid->nodes[i].coords[0];
            p[1]=grid->nodes[i].coords[1];
            p[2]=grid->nodes[i].coords[2];
            ignore=fwrite(p, sizeof(double), 3, F);
        }
        for(i=0; i<grid->NN; i++)
        {
            ignore=fwrite(&grid->nodes[i].on_bdry, sizeof(int), 1, F);
        }
    }
    ignore=fwrite(&grid->NEL, sizeof(int), 1, F);
    {
        int i;
        for(i=0; i<grid->NEL; i++)
        {
            ignore=fwrite(grid->elems[i].nodes, sizeof(int), 4, F);
        }
    }
    if(ferror(F))
    {
        perror("SaveGridBin");
        error_msg("SaveGridBin: fwrite failed.\n");
    }
    fclose(F);
    return;
}


/****************************************************************************/
/*            LoadGrid */

/****************************************************************************/

void MakeGridNameFromPath(grid3d_t *grid, const char *path)
{
    char const *fname;
    char *s;
    int len;
    fname=strrchr(path, '/');
    if(fname==NULL)
        fname=path;
    else
        fname++;
    len=strlen(fname);
    grid->gridname=calloc(len+1, sizeof(char));
    strncpy(grid->gridname, fname, len);
    if((s=strrchr(grid->gridname, '.'))!=NULL)
        *s=0;
}

grid3d_t *LoadGridBin(const char *path)
{
    FILE *F;
    grid3d_t *grid;
    int i;

    LOG_MPI_WTIME_START;
    F=fopen(path, "rb");
    if(F==NULL)
    {
        perror("LoadGridBin");
        error_msg("LoadGridBin: cannot open file %s for reading.\n", path);
    }

    grid=(grid3d_t *)calloc(1, sizeof(grid3d_t));
    if(grid==NULL)
        goto alloc_error;

    /* read nodes */
    if(fread(&grid->NN, sizeof(int), 1, F)<1)
        goto read_error;
    grid->nodes=malloc(grid->NN*sizeof(node_t));
    if(grid->nodes==NULL)
        goto alloc_error;
    for(i=0; i<grid->NN; i++)
    {
        double X[3];
        node_t *n=grid->nodes+i;
        if(fread(X, sizeof(double), 3, F)<3)
            goto read_error;
        n->coords[0]=X[0];
        n->coords[1]=X[1];
        n->coords[2]=X[2];
    }
    grid->NNI=0;
    for(i=0; i<grid->NN; i++)
    {
        node_t *n=grid->nodes+i;
        if(fread(&n->on_bdry, sizeof(unsigned), 1, F)<1)
            goto read_error;
        if(n->on_bdry)
            n->bcond=BC_DIR;
        else
            grid->NNI++;
    }
    /* read elements */
    if(fread(&grid->NEL, sizeof(int), 1, F)<1)
        goto read_error;
    grid->elems=malloc(grid->NEL*sizeof(elem_t));
    if(grid->nodes==NULL)
        goto alloc_error;
    for(i=0; i<grid->NEL; i++)
    {
        if(fread(grid->elems[i].nodes, sizeof(int), 4, F)<4)
            goto read_error;
    }

    fclose(F);

    EnrichGrid(grid);
    UpdateMeasures(grid);
    MakeGridNameFromPath(grid, path);

    LOG_MPI_WTIME_STOP;
    LogGridStats(grid);
    return grid;

alloc_error:
    error_msg("LoadGridBin: memory allocation error.\n");
read_error:
    error_msg("LoadGridBin: cannot read from file %s.\n", path);
}

grid3d_t *LoadGrid(const char *path)
{
    FILE *F;
    grid3d_t *grid;
    int i, ignor;

    LOG_MPI_WTIME_START;
    F=fopen(path, "rt");
    if(F==NULL)
    {
        perror("LoadGrid");
        error_msg("LoadGrid: cannot open file %s\n", path);
    }

    grid=calloc(1, sizeof(grid3d_t));
    if(grid==NULL)
        goto alloc_error;

    ignor=fscanf(F, "%u", &grid->NN);
    while(fgetc(F)!='\n');

    grid->nodes=calloc(grid->NN, sizeof(node_t));
    if(grid->nodes==NULL)
        goto alloc_error;

    for(i=0; i<grid->NN; i++)
    {
        double p[3];
        node_t *n=grid->nodes+i;
        ignor=fscanf(F, "%lg %lg %lg", p, p+1, p+2);
        n->coords[0]=p[0];
        n->coords[1]=p[1];
        n->coords[2]=p[2];
    }

    grid->NNI=0;
    for(i=0; i<grid->NN; i++)
    {
        node_t *n=grid->nodes+i;
        ignor=fscanf(F, "%u", &n->on_bdry);
        if(n->on_bdry)
            n->bcond=BC_DIR;
        else
            grid->NNI++;
    }

    ignor=fscanf(F, "%u", &grid->NEL);
    while(fgetc(F)!='\n');

    grid->elems=calloc(grid->NEL, sizeof(elem_t));
    if(grid->elems==NULL)
        goto alloc_error;

    for(i=0; i<grid->NEL; i++)
    {
        elem_t *T=grid->elems+i;
        ignor=fscanf(F, "%u %u %u %u", T->nodes, T->nodes+1, T->nodes+2, T->nodes+3);
    }

    fclose(F);

    EnrichGrid(grid);
    UpdateMeasures(grid);
    MakeGridNameFromPath(grid, path);

    LOG_MPI_WTIME_STOP;
    LogGridStats(grid);
    return grid;

alloc_error:
    error_msg("LoadGrid: memory allocation error.\n");
}

/****************************************************************************/

void save_grid_bin(const grid3d_t *grid, FILE *F)
{
    if(mpi_my_id!=0)
        return;

    int i; size_t ignore;
    LOG_MPI_WTIME_START;

    if(F==NULL) return;
    ignore=fwrite(&grid->NN, sizeof(int), 1, F);
    ignore=fwrite(&grid->NNI, sizeof(int), 1, F);
    ignore=fwrite(&grid->NF, sizeof(int), 1, F);
    ignore=fwrite(&grid->NFI, sizeof(int), 1, F);
    ignore=fwrite(&grid->NE, sizeof(int), 1, F);
    ignore=fwrite(&grid->NEI, sizeof(int), 1, F);
    ignore=fwrite(&grid->NDFI, sizeof(int), 1, F);
    ignore=fwrite(&grid->NEL, sizeof(int), 1, F);
    for(i=0; i<grid->NN; i++)
    {
        node_t *nd=grid->nodes+i;
        ignore=fwrite(nd->coords, sizeof(double), 3, F);
        ignore=fwrite(&nd->bcond, sizeof(int), 1, F);
        ignore=fwrite(&nd->on_bdry, sizeof(int), 1, F);
    }
    for(i=0; i<grid->NE; i++)
    {
        edge_t *ed=grid->edges+i;
        ignore=fwrite(ed->nodes, sizeof(int), 2, F);
    }
    for(i=0; i<grid->NF; i++)
    {
        face_t *fc=grid->faces+i;
        ignore=fwrite(fc->nodes, sizeof(int), 3, F);
        ignore=fwrite(fc->edges, sizeof(int), 3, F);
        ignore=fwrite(fc->elems, sizeof(int), 2, F);
    }
    for(i=0; i<grid->NEL; i++)
    {
        elem_t *el=grid->elems+i;
        ignore=fwrite(el->nodes, sizeof(int), 4, F);
        ignore=fwrite(el->edges, sizeof(int), 6, F);
        ignore=fwrite(el->faces, sizeof(int), 4, F);
    }
    if(ferror(F)) perror(NULL);
    LOG_MPI_WTIME_STOP;
    return;
}

void load_grid_bin(grid3d_t *grid, FILE *F)
{
    int i;
    LOG_MPI_WTIME_START;

    if(F==NULL) return;
    if(fread(&grid->NN, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NNI, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NF, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NFI, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NE, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NEI, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NDFI, sizeof(int), 1, F)<1)
        goto read_error;
    if(fread(&grid->NEL, sizeof(int), 1, F)<1)
        goto read_error;
    grid->nodes=malloc(grid->NN*sizeof(node_t));
    if(grid->nodes==NULL)
        goto alloc_error;
    for(i=0; i<grid->NN; i++)
    {
        node_t *nd=grid->nodes+i;
        if(fread(nd->coords, sizeof(double), 3, F)<3)
            goto read_error;
        if(fread(&nd->bcond, sizeof(unsigned), 1, F)<1)
            goto read_error;
        if(fread(&nd->on_bdry, sizeof(unsigned), 1, F)<1)
            goto read_error;
    }
    grid->edges=malloc(grid->NE*sizeof(edge_t));
    if(grid->edges==NULL)
        goto alloc_error;
    for(i=0; i<grid->NE; i++)
    {
        edge_t *ed=grid->edges+i;
        if(fread(ed->nodes, sizeof(int), 2, F)<2)
            goto read_error;
    }
    grid->faces=malloc(grid->NF*sizeof(face_t));
    if(grid->faces==NULL)
        goto alloc_error;
    for(i=0; i<grid->NF; i++)
    {
        face_t *fc=grid->faces+i;
        if(fread(fc->nodes, sizeof(int), 3, F)<3)
            goto read_error;
        if(fread(fc->edges, sizeof(int), 3, F)<3)
            goto read_error;
        if(fread(fc->elems, sizeof(int), 2, F)<2)
            goto read_error;
    }
    grid->elems=malloc(grid->NEL*sizeof(elem_t));
    if(grid->elems==NULL)
        goto alloc_error;
    for(i=0; i<grid->NEL; i++)
    {
        elem_t *el=grid->elems+i;
        if(fread(el->nodes, sizeof(int), 4, F)<4)
            goto read_error;
        if(fread(el->edges, sizeof(int), 6, F)<6)
            goto read_error;
        if(fread(el->faces, sizeof(int), 4, F)<4)
            goto read_error;
    }
    if(ferror(F)) perror(NULL);

    LOG_MPI_WTIME_STOP;
    LogGridStats(grid);

    return;

alloc_error:
    error_msg("%s: memory allocation error.\n", __func__);
read_error:
    error_msg("%s: cannot read from file.\n", __func__);
}



